/*!
 * \file gico_object.c
 *
 */

#include "gico_object.h"

/* TODO Check the declared variables of each function and if they are being used or not */

int gico_object_alloc(gico_object **newone, size_t N1, size_t N2, enum gico_object_kind kind)
{
	/* TODO: Check that the object is not previously initialized to avoid memory leaks */
	int ret = GSL_SUCCESS;		/* returning value */
	size_t size;				/* size of the data */

	/* allocation of the gico_object */
	*newone = calloc(1, sizeof(gico_object));
	if (*newone == NULL)
	{
		ret = GSL_ENOMEM;
		goto out_nomemory;
	}

	/* allocation of the data vector */
	if (kind == GICO_COMPLEX)
	{
		size = (size_t) 2*N1*N2;
	}
	else
	{
		size = (size_t) N1*N2;
	}
	(*newone)->data = calloc(size, sizeof(double));

	if ((*newone)->data == NULL)
	{
		ret = GSL_ENOMEM;
		goto out_nomemory;
	}

	/* initializing the variables */
	(*newone)->N[0] = N1;
	(*newone)->N[1] = N2;
	(*newone)->kind = kind;

    ret = GSL_SUCCESS;
	goto out;

out_nomemory:

	if (newone != NULL)
	{
		free(newone);
	}

out:

	return ret;
}


int gico_object_dealloc(gico_object *oldone)
{
	int ret = GSL_SUCCESS;		/* returning value */

	if (oldone->data != NULL)
	{
		free(oldone->data);
	}

	if (oldone != NULL)
	{
		free(oldone);
	}

	ret = GSL_SUCCESS;
	goto out;

out:

	return ret;
}


int gico_object_set_value(gico_object *obj, size_t n1, size_t n2, double value)
{
	int ret = GSL_SUCCESS;		/* returning value */

	/* check if the object is initialized */
	if (!GICO_OBJECT_CHECK_INIT(obj))
	{
		ret = GSL_EINVAL;
		goto out_noobj;
	}

	/* check if the position is inside the boundaries */
	if ((n1 >= obj->N[0]) || (n1 < 0) || (n2 >= obj->N[1]) || (n2 < 0))
	{
		ret = GSL_ERANGE;
		goto out_outbounds;
	}

	/* insert the value */
	obj->data[(size_t) n1*obj->N[0]+n2] = value;

	ret = GSL_SUCCESS;
	goto out;

out_outbounds:

out_noobj:

out:

	return ret;
}


int gico_object_get_value(gico_object *obj, size_t n1, size_t n2, double *value)
{
	int ret = GSL_SUCCESS;		/* returning value */

	/* check that the object is initialized */
	if (!GICO_OBJECT_CHECK_INIT(obj))
	{
		ret = GSL_EINVAL;
		goto out_noobj;
	}

	/* check that the value is in range */
	if ((n1 >= obj->N[0]) || (n1 < 0) || (n2 >= obj->N[1]) || (n2 < 0))
	{
		ret = GSL_ERANGE;
		goto out_outbounds;
	}

	/* get the value */
	*value = obj->data[(size_t) n1*obj->N[0]+n2];

	ret = GSL_SUCCESS;
	goto out;

out_outbounds:

out_noobj:

out:

	return ret;
}


int gico_object_set_cvalue(gico_object *obj, size_t n1, size_t n2, gsl_complex cvalue)
{
	int ret = GSL_SUCCESS;		/* returning value */

	/* check if the object is initialized */
	if (!GICO_OBJECT_CHECK_INIT(obj))
	{
		ret = GSL_EINVAL;
		goto out_noobj;
	}

	/* TODO: Check that obj is complex? */

	/* check if the position is inside the boundaries */
	if ((n1 >= obj->N[0]) || (n1 < 0) || (n2 >= obj->N[1]) || (n2 < 0))
	{
		ret = GSL_ERANGE;
		goto out_outbounds;
	}

	/* insert the value */
	obj->data[(size_t) 2*(n1*obj->N[0]+n2)] = cvalue.dat[0];
	obj->data[(size_t) 2*(n1*obj->N[0]+n2)+1] = cvalue.dat[1];

	ret = GSL_SUCCESS;
	goto out;

out_outbounds:

out_noobj:

out:

	return ret;
}


int gico_object_get_cvalue(gico_object *obj, size_t n1, size_t n2, gsl_complex *cvalue)
{
	int ret = GSL_SUCCESS;		/* returning value */

	/* check that the object is initialized */
	if (!GICO_OBJECT_CHECK_INIT(obj))
	{
		goto out_noobj;
	}

	/* TODO: Check that obj is complex? */

	/* check that the value is in range */
	if ((n1 >= obj->N[0]) || (n1 < 0) || (n2 >= obj->N[1]) || (n2 < 0))
	{
		goto out_outbounds;
	}

	/* get the value */
	cvalue->dat[0] = obj->data[(size_t) 2*(n1*obj->N[0]+n2)];
	cvalue->dat[1] = obj->data[(size_t) 2*(n1*obj->N[0]+n2)+1];

	ret = GSL_SUCCESS;
	goto out;

out_outbounds:

	if (ret == GSL_SUCCESS)
	{
		ret = GSL_ERANGE;
	}

out_noobj:

	if (ret == GSL_SUCCESS)
	{
		ret = GSL_EINVAL;
	}

out:

	return ret;
}


int gico_object_set_vector(gico_object *obj, size_t n1, size_t n2, enum gico_vector_align va, gsl_vector *vector)
{
	int ret = GSL_SUCCESS;		/* returning value */
	int i;						/* iteration token */

	/* check if the object is initialized */
	if (!GICO_OBJECT_CHECK_INIT(obj))
	{
		ret = GSL_EINVAL;
		goto out_noobj;
	}

	/* check if the position is inside the boundaries */
	if (((va == GICO_VER) && ((n1 >= obj->N[0] - vector->size) || (n1 < 0) || (n2 >= obj->N[1]) || (n2 < 0))) ||
		((va == GICO_HOR) && ((n1 >= obj->N[0]) || (n1 < 0) || (n2 >= obj->N[1] - vector->size) || (n2 < 0))))
	{
		ret = GSL_ERANGE;
		goto out_outbounds;
	}

	/* insert the values of the vector */
	if (va == GICO_VER)
	{
		for (i=0; i<vector->size; i++)
		{
			obj->data[(size_t) (n1+i)*obj->N[0]+n2] = vector->data[i];
		}
	}
	else
	{
		for (i=0; i<vector->size; i++)
		{
			obj->data[(size_t) n1*obj->N[0]+n2+i] = vector->data[i];
		}
	}

	ret = GSL_SUCCESS;
	goto out;

out_outbounds:

out_noobj:

out:

	return ret;
}


int gico_object_get_vector(gico_object *obj, size_t n1, size_t n2, enum gico_vector_align va, gsl_vector *vector)
{
	int ret = GSL_SUCCESS;		/* returning value */
	int i;						/* iteration token */

	/* check if the object is initialized */
	if (!GICO_OBJECT_CHECK_INIT(obj))
	{
		ret = GSL_EINVAL;
		goto out_noobj;
	}

	/* check if the position is inside the boundaries */
	if (((va == GICO_VER) && ((n1 >= obj->N[0] - vector->size) || (n1 < 0) || (n2 >= obj->N[1]) || (n2 < 0))) ||
		((va == GICO_HOR) && ((n1 >= obj->N[0]) || (n1 < 0) || (n2 >= obj->N[1] - vector->size) || (n2 < 0))))
	{
		ret = GSL_ERANGE;
		goto out_outbounds;
	}

	/* get the values for the vector */
	if (va == GICO_VER)
	{
		for (i=0; i<vector->size; i++)
		{
			vector->data[i] = obj->data[(size_t) (n1+i)*obj->N[0]+n2];
		}
	}
	else
	{
		for (i=0; i<vector->size; i++)
		{
			 vector->data[i] = obj->data[(size_t) n1*obj->N[0]+n2+i];
		}
	}

	ret = GSL_SUCCESS;
	goto out;

out_outbounds:

out_noobj:

out:

	return ret;
}


int gico_object_set_cvector(gico_object *obj, size_t n1, size_t n2, enum gico_vector_align va, gsl_vector_complex *vector)
{
	int i=0;
	int ret = GSL_SUCCESS;		/* returning value */

	/* check if the object is initialized */
	if (!GICO_OBJECT_CHECK_INIT(obj))
	{
		ret = GSL_EINVAL;
		goto out_noobj;
	}

	/* check if the position is inside the boundaries */
	if (((va == GICO_VER) && ((n1 >= obj->N[0] - vector->size) || (n1 < 0) || (n2 >= obj->N[1]) || (n2 < 0))) ||
		((va == GICO_HOR) && ((n1 >= obj->N[0]) || (n1 < 0) || (n2 >= obj->N[1] - vector->size) || (n2 < 0))))
	{
		ret = GSL_ERANGE;
		goto out_outbounds;
	}

	/* insert the values of the vector */
	if (va == GICO_VER)
	{
		for (i=0; i<vector->size; i++)
		{
			obj->data[(size_t) 2*((n1+i)*obj->N[0]+n2)] = vector->data[(size_t) 2*i];
			obj->data[(size_t) 2*((n1+i)*obj->N[0]+n2)+1] = vector->data[(size_t) 2*i+1];
		}
	}
	else
	{
		for (i=0; i<vector->size; i++)
		{
			obj->data[(size_t) 2*(n1*obj->N[0]+n2+i)] = vector->data[(size_t) 2*i];
			obj->data[(size_t) 2*(n1*obj->N[0]+n2+i)+1] = vector->data[(size_t) 2*i+1];
		}
	}

	ret = GSL_SUCCESS;
	goto out;

out_outbounds:

out_noobj:

out:

	return ret;
}


int gico_object_get_cvector(gico_object *obj, size_t n1, size_t n2, enum gico_vector_align va, gsl_vector_complex *vector)
{
	int i=0;
	int ret = GSL_SUCCESS;		/* returning value */

	/* check if the object is initialized */
	if (!GICO_OBJECT_CHECK_INIT(obj))
	{
		ret = GSL_EINVAL;
		goto out_noobj;
	}

	/* check if the position is inside the boundaries */
	if (((va == GICO_VER) && ((n1 >= obj->N[0] - vector->size) || (n1 < 0) || (n2 >= obj->N[1]) || (n2 < 0))) ||
		((va == GICO_HOR) && ((n1 >= obj->N[0]) || (n1 < 0) || (n2 >= obj->N[1] - vector->size) || (n2 < 0))))
	{
		ret = GSL_ERANGE;
		goto out_outbounds;
	}

	/* insert the values of the vector */
	if (va == GICO_VER)
	{
		for (i=0; i<vector->size; i++)
		{
			vector->data[(size_t) 2*i] = obj->data[(size_t) 2*((n1+i)*obj->N[0]+n2)];
			vector->data[(size_t) 2*i+1] = obj->data[(size_t) 2*((n1+i)*obj->N[0]+n2)+1];
		}
	}
	else
	{
		for (i=0; i<vector->size; i++)
		{
			vector->data[(size_t) 2*i] = obj->data[(size_t) 2*(n1*obj->N[0]+n2+i)];
			vector->data[(size_t) 2*i+1] = obj->data[(size_t) 2*(n1*obj->N[0]+n2+i)+1];
		}
	}

	ret = GSL_SUCCESS;
	goto out;

out_outbounds:

out_noobj:

out:

	return ret;
}


int gico_object_set_matrix(gico_object *obj, size_t n1, size_t n2, gsl_matrix *matrix)
{
	int i=0, j=0;
	size_t mtda = 0;
	int ret = GSL_SUCCESS;		/* returning value */

	/* check if the object is initialized */
	if (!GICO_OBJECT_CHECK_INIT(obj))
	{
		ret = GSL_EINVAL;
		goto out_noobj;
	}

	/* check if the position is inside the boundaries */
	if ((n1 >= obj->N[0] - matrix->size1) || (n1 < 0) || (n2 >= obj->N[1] - matrix->size2) || (n2 < 0))
	{
		ret = GSL_ERANGE;
		goto out_outbounds;
	}

	/* get the matrix tda */
	mtda = matrix->tda;

	/* insert the values of the vector */
	for (i=0; i<matrix->size1; i++)
	{
		for (j=0; j<matrix->size2; j++)
		{
			obj->data[(size_t) (n1+i)*obj->N[0]+n2+j] = matrix->data[(size_t) i*mtda+j];
		}
	}

	ret = GSL_SUCCESS;
	goto out;

out_outbounds:

out_noobj:

out:

	return ret;
}


int gico_object_get_matrix(gico_object *obj, size_t n1, size_t n2, gsl_matrix *matrix)
{
	int i=0, j=0;
	size_t mtda = 0;
	int ret = GSL_SUCCESS;		/* returning value */

	/* check if the object is initialized */
	if (!GICO_OBJECT_CHECK_INIT(obj))
	{
		ret = GSL_EINVAL;
		goto out_noobj;
	}

	/* check if the position is inside the boundaries */
	if ((n1 >= obj->N[0] - matrix->size1) || (n1 < 0) || (n2 >= obj->N[1] - matrix->size2) || (n2 < 0))
	{
		ret = GSL_ERANGE;
		goto out_outbounds;
	}

	/* get the matrix tda */
	mtda = matrix->tda;

	/* insert the values of the vector */
	for (i=0; i<matrix->size1; i++)
	{
		for (j=0; j<matrix->size2; j++)
		{
			matrix->data[(size_t) i*mtda+j] = obj->data[(size_t) (n1+i)*obj->N[0]+n2+j];
		}
	}

	ret = GSL_SUCCESS;
	goto out;

out_outbounds:

out_noobj:

out:

	return ret;
}


int gico_object_set_cmatrix(gico_object *obj, size_t n1, size_t n2, gsl_matrix_complex *matrix)
{
	int i=0, j=0;
	size_t mtda = 0;
	int ret = GSL_SUCCESS;		/* returning value */

	/* check if the object is initialized */
	if (!GICO_OBJECT_CHECK_INIT(obj))
	{
		ret = GSL_EINVAL;
		goto out_noobj;
	}

	/* check if the position is inside the boundaries */
	if ((n1 >= obj->N[0] - matrix->size1) || (n1 < 0) || (n2 >= obj->N[1] - matrix->size2) || (n2 < 0))
	{
		ret = GSL_ERANGE;
		goto out_outbounds;
	}

	/* get the matrix tda */
	mtda = matrix->tda;

	/* insert the values of the vector */
	for (i=0; i<matrix->size1; i++)
	{
		for (j=0; j<matrix->size2; j++)
		{
			obj->data[(size_t) 2*((n1+i)*obj->N[0]+n2+j)] = matrix->data[(size_t) 2*(i*mtda+j)];
			obj->data[(size_t) 2*((n1+i)*obj->N[0]+n2+j)+1] = matrix->data[(size_t) 2*(i*mtda+j)+1];
		}
	}

	ret = GSL_SUCCESS;
	goto out;

out_outbounds:

out_noobj:

out:

	return ret;
}


int gico_object_get_cmatrix(gico_object *obj, size_t n1, size_t n2, gsl_matrix_complex *matrix)
{
	int i=0, j=0;
	size_t mtda = 0;
	int ret = GSL_SUCCESS;		/* returning value */

	/* check if the object is initialized */
	if (!GICO_OBJECT_CHECK_INIT(obj))
	{
		ret = GSL_EINVAL;
		goto out_noobj;
	}

	/* check if the position is inside the boundaries */
	if ((n1 >= obj->N[0] - matrix->size1) || (n1 < 0) || (n2 >= obj->N[1] - matrix->size2) || (n2 < 0))
	{
		ret = GSL_ERANGE;
		goto out_outbounds;
	}

	/* get the matrix tda */
	mtda = matrix->tda;

	/* insert the values of the vector */
	for (i=0; i<matrix->size1; i++)
	{
		for (j=0; j<matrix->size2; j++)
		{
			matrix->data[(size_t) 2*(i*mtda+j)] = obj->data[(size_t) 2*((n1+i)*obj->N[0]+n2+j)];
			matrix->data[(size_t) 2*(i*mtda+j)+1] = obj->data[(size_t) 2*((n1+i)*obj->N[0]+n2+j)+1];
		}
	}

	ret = GSL_SUCCESS;
	goto out;

out_outbounds:

out_noobj:

out:

	return ret;
}


int gico_object_save(gico_object *obj, const char filename[])
{
	int factor = 1;				/* multiplicative factor in case of complex data */
	int fd;						/* file descriptor */
	int ret = GSL_SUCCESS;		/* returning value */
	int i,j;
	size_t n_bytes, n_write;

	/* check if the object is initialized */
	if (!GICO_OBJECT_CHECK_INIT(obj))
	{
		ret = GSL_EINVAL;
		goto out_noobj;
	}

    /* if the file exists delete it */
	/* TODO: Implement an asking version? */
    if (access(filename, F_OK) != -1)
    {
        /* delete the file */
        ret = unlink(filename);

        /* if unable to delete the file */
        if(ret == -1)
        {
        	ret = GSL_EFAILED;
        	goto out_nodelete;
        }
    }

    /* open the file for new creation */
    fd = open(filename, O_RDWR | O_APPEND | O_CREAT,
                  S_IRUSR | S_IRGRP | S_IWUSR | S_IWGRP | S_IROTH);

    /* check if we were able to open the file */
    if (!fd)
    {
    	ret = GSL_EFAULT;
        goto out_noopen;
    }

    /* dump the gico_signal dimensions */
    ret = gico_write(fd, (const size_t *) (obj->N), (size_t) 2*sizeof(size_t));
    if (ret != GSL_SUCCESS)
    {
    	ret = GSL_EOF;
    	goto out_nowrite;
    }

    /* dump the physical dimensions of the matrix */
    ret = gico_write(fd, (const double *) (obj->d), (size_t) 2*sizeof(double));
    if (ret != GSL_SUCCESS)
    {
    	ret = GSL_EOF;
    	goto out_nowrite;
    }

    /* dump the kind of object */
    ret = gico_write(fd, (const int *) &(obj->kind), sizeof(int));
    if (ret != GSL_SUCCESS)
    {
    	ret = GSL_EOF;
    	goto out_nowrite;
    }

    /* get the right number of doubles per data */
    if (obj->kind == GICO_COMPLEX)
    {
    	factor = 2;
    }

    /* dump the data */
    ret = gico_write(fd, (const char *) (obj->data), (size_t) factor*obj->N[0]*obj->N[1]*sizeof(double));
    if (ret != GSL_SUCCESS)
    {
    	ret = GSL_EOF;
    	goto out_nowrite;
    }

    /* Close the file */
    close(fd);

    ret = GSL_SUCCESS;
    goto out;

out_nowrite:
	close(fd);

out_noopen:

out_nodelete:

out_noobj:

out:

	return ret;
}

int gico_object_load(const char filename[], gico_object **obj)
{
	int factor = 1;				/* multiplicative factor in case of complex data */
	int fd;						/* file descriptor */
	int ret = GSL_SUCCESS;		/* returning value */
	int i, j;
	size_t n_bytes, n_read;

	/* temporal variables to store the data */
	size_t N[2];
	double d[2];
	enum gico_object_kind kind;

    /* if the file do not exists exit */
    if (access(filename, F_OK) == -1)
    {
    	ret = GSL_EFAILED;
        goto out_noexists;
    }

    /* open the file for reading */
    fd = open(filename, O_RDONLY, S_IRUSR | S_IRGRP | S_IWUSR | S_IWGRP | S_IROTH);

    /* check if we were able to open the file */
    if (fd < 0)
    {
    	ret = GSL_EFAULT;
        goto out_noopen;
    }

    /* dump the gico_signal dimensions */
    ret = gico_read(fd, (size_t *) N, (size_t) 2*sizeof(size_t));
    if (ret != GSL_SUCCESS)
    {
    	ret = GSL_EOF;
    	goto out_nowrite;
    }

    /* dump the physical dimensions of the matrix */
    ret = gico_read(fd, (double *) d, (size_t) 2*sizeof(double));
    if (ret != GSL_SUCCESS)
    {
    	ret = GSL_EOF;
    	goto out_nowrite;
    }

    /* dump the kind of object */
    ret = gico_read(fd, (int *) &kind, sizeof(int));
    if (ret != GSL_SUCCESS)
    {
    	ret = GSL_EOF;
    	goto out_nowrite;
    }

    /* initialize the gico_object */
    ret = gico_object_alloc(obj, N[0], N[1], kind);
    if (ret != GSL_SUCCESS)
    {
    	ret = GSL_EOF;
    	goto out_noobjalloc;
    }

    /* insert the known values */
    (*obj)->d[0] = d[0];
    (*obj)->d[1] = d[1];

    /* get the right number of doubles per data */
    if ((*obj)->kind == GICO_COMPLEX)
    {
    	factor = 2;
    }

    /* dump the data */
    ret = gico_read(fd, (double *) ((*obj)->data), (size_t) factor*(*obj)->N[0]*(*obj)->N[1]*sizeof(double));
    if (ret != GSL_SUCCESS)
    {
    	ret = GSL_EOF;
    	goto out_nowrite;
    }

    /* close the file */
    close(fd);

    ret = GSL_SUCCESS;
    goto out;

out_noobjalloc:
out_nowrite:
	close(fd);

out_noopen:

out_noexists:

out_noobj:

out:

	return ret;
}


int gico_object_save_gnuplot(gico_object *obj, enum gico_object_kind kind, const char filename[])
{
	int factor = 1;				/* multiplicative factor in case of complex data */
	int fd;						/* file descriptor */
	int ret = GSL_SUCCESS;		/* returning value */
	int i,j;
	float Nx, Ny;
	float d;
	gsl_complex c;

	/* check if the object is initialized */
	if (!GICO_OBJECT_CHECK_INIT(obj))
	{
		ret = GSL_EINVAL;
		goto out_noobj;
	}

    /* if the file exists delete it */
	/* TODO: Implement an asking version? */
    if (access(filename, F_OK) != -1)
    {
        /* delete the file */
        ret = unlink(filename);

        /* if unable to delete the file */
        if(ret == -1)
        {
        	ret = GSL_EFAILED;
        	goto out_nodelete;
        }
    }

    /* open the file for new creation */
    fd = open(filename, O_RDWR | O_APPEND | O_CREAT,
                  S_IRUSR | S_IRGRP | S_IWUSR | S_IWGRP | S_IROTH);

    /* check if we were able to open the file */
    if (!fd)
    {
    	ret = GSL_EFAULT;
        goto out_noopen;
    }

    /* save the first row */
    Nx = (float) (obj->N[1]);
    ret = gico_write(fd, (const float *) &Nx, (size_t) sizeof(float));
    if (ret != GSL_SUCCESS)
    {
    	ret = GSL_EOF;
    	goto out_nowrite;
    }

    for (i=0; i<Nx; i++)
    {
		d = (float) (-ceil(Nx/2)+1+i)*obj->d[1];
	    ret = gico_write(fd, (const float *) &d, sizeof(float));
	    if (ret != GSL_SUCCESS)
	    {
	    	ret = GSL_EOF;
	    	goto out_nowrite;
	    }
    }


    /* save each row */
    Ny = (double) (obj->N[0]);
    for (i=0; i<Ny; i++)
    {
		d = (float) (-ceil(Ny/2)+1+i)*obj->d[0];
	    ret = gico_write(fd, (const float *) &d, sizeof(float));
	    if (ret != GSL_SUCCESS)
	    {
	    	ret = GSL_EOF;
	    	goto out_nowrite;
	    }

	    for (j=0; j<Nx; j++)
	    {
	    	if (obj->kind == GICO_COMPLEX)
	    	{
	    		switch(kind)
				{
				case GICO_REAL:
					d = (float) obj->data[2*(j + i*(size_t)Nx)];
					break;
				case GICO_IMAG:
					d = (float) obj->data[1 + 2*(j + i*(size_t)Nx)];
					break;
				case GICO_ABS:
				case GICO_ABS2:
				case GICO_ANGLE:
					GSL_SET_REAL(&c, obj->data[2*(j + i*(size_t)Nx)]);
					GSL_SET_IMAG(&c, obj->data[1 + 2*(j + i*(size_t)Nx)]);
					if (kind == GICO_ABS)
					{
						d = (float) gsl_complex_abs(c);
					}
					else if (kind == GICO_ABS2)
					{
						d = (float) gsl_complex_abs2(c);
					}
					else
					{
						d = (float) gsl_complex_arg(c);
					}
					break;

				default:
					ret = GSL_EINVAL;
					goto out_nowrite;
				}
	    	}
			else
			{
				d = (float) obj->data[j + i*(size_t)Nx];
				switch(kind)
				{
				case GICO_REAL:
					break;
				case GICO_IMAG:
					d = (float) 0.0;
					break;
				case GICO_ABS:
					d = (float) fabsf(d);
					break;
				case GICO_ABS2:
					d = (float) pow(fabs(d), 2.0);
					break;
				case GICO_ANGLE:
					GSL_SET_COMPLEX(&c, d, 0.0);
					d = (float) gsl_complex_arg(c);
					break;

				default:
					ret = GSL_EINVAL;
					goto out_nowrite;
				}
			}

	    	ret = gico_write(fd, (const float *) &d, sizeof(float));
	    	if (ret != GSL_SUCCESS)
			{
				ret = GSL_EOF;
				goto out_nowrite;
			}
	    }
    }

    /* make sure we have written everything */
    fsync(fd);

    /* close the file */
    close(fd);

    ret = GSL_SUCCESS;
    goto out;

out_nowrite:
	close(fd);

out_noopen:

out_nodelete:

out_complex:

out_noobj:

out:

	return ret;
}


int gico_write(int fd, const void *buffer, size_t n_bytes)
{
    int ret = 0;
    ssize_t n_write = 0;

    do
    {
    	n_write = write(fd, (const char *) (buffer+ret), n_bytes-ret);
    	ret += n_write;
    } while ((ret != n_bytes) && (n_write > 0));

    if (n_write < 0)
    {
    	ret = GSL_FAILURE;
    }
    else
    {
    	ret = GSL_SUCCESS;
    }

    return ret;
}


int gico_read(int fd, void *buffer, size_t n_bytes)
{
    size_t ret;
    size_t n_read;

    n_read = read(fd, (char *) (buffer+ret), n_bytes-ret);

    if (n_read < n_bytes)
    {
    	ret = GSL_FAILURE;
    }
    else
    {
    	ret = GSL_SUCCESS;
    }

    return ret;
}
